arXiv: 1501.02912v2 [cond-mat.soft] 17 Dec 2015 


Large-scale simulation of steady and time-dependent active 
suspensions with the force-coupling method 


Blaise Delmotte'^’'’’*, Eric E. Keaveny'^’*, Franck Plouraboue’^’'^, Eric Climent'^’’’ 

“ University of Toulouse INPT- UPS: Institut de Mecanique des Fluides, Toulouse, France. 

^IMFT - CNRS, UMR 5502 1, Alice du Professeur Camille Soula, 31 400 Toulouse, France. 
‘^Department of Mathematics, Imperial College London, South Kensington Campus, London, SW7 2AZ, 

UK 


Abstract 

We present a new development of the force-conpling method (FCM) to address the accn- 
rate simnlation of a large nnmber of interacting micro-swimmers. Onr approach is based on 
the sqnirmer model, which we adapt to the FCM framework, resnlting in a method that is 
snitable for simnlating semi-dilnte sqnirmer snspensions. Other effects, snch as steric inter¬ 
actions, are considered with onr model. We test onr method by comparing the velocity held 
aronnd a single sqnirmer and the pairwise interactions between two sqnirmers with exact 
solntions to the Stokes eqnations and results given by other numerical methods. We also il¬ 
lustrate our method’s ability to describe spheroidal swimmer shapes and biologically-relevant 
time-dependent swimming gaits. We detail the numerical algorithm used to compute the 
hydrodynamic coupling between a large collection (10^ — 10^) of micro-swimmers. Using this 
methodology, we investigate the emergence of polar order in a suspension of squirmers and 
show that for large domains, both the steady-state polar order parameter and the growth 
rate of instability are independent of system size. These results demonstrate the effective¬ 
ness of our approach to achieve near continuum-level results, allowing for better comparison 
with experimental measurements while complementing and informing continuum models. 

Keywords: Force Coupling Method, low Reynolds number, active suspension, swimming 
gait, collective dynamics. High Performance Computing 


1. Introduction 

Suspensions of active, self-propelled particles arise in both biological systems, such as 
populations of micro-organisms [DEI El ID and synthetic, colloidal systems [S]. These sus¬ 
pensions can exhibit the formation of coherent structures and complex flow patterns which 
may lead to enhanced mixing of chemicals in the surrounding fluid, the alteration of sus¬ 
pension rheology, or, in the biological case, increased nutrient uptake by a population of 
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micro-organisms. In addition to promising applications such as algae biofuels EE], charac¬ 
terizing the collective dynamics found in these suspensions is of fundamental importance to 
understanding zooplankton dynamics and mammal fertility mm- 

The mathematical modeling of active suspensions entails describing how individual swim¬ 
mers move and interact in response to the flow helds that they generate [niEallll]. It is 
particularly important for these models to be able to handle a large collection of swim¬ 
mers in order to obtain suspension properties at the lab/m situ scale. The modeling of 
the collective behavior of active matter has been a vibrant area of research during the last 
decade [laiiHiiiiiizi, to cite only a few recent reviews. Generally speaking, the modelling 
approaches can be sorted into two categories: continuum theories and particle-based simu¬ 
lations. Most of the continuum models are generally valid for dilute suspensions where the 
hydrodynamic disturbances are given by a mean-held description of far-held hydrodynamic 
interactions [HI Eg EE]. Recent advances towards more concentrated suspensions include 
steric interactions IZDl, but the inclusion of high-order singularities due to particle size re¬ 
mains outstanding. Despite this, these models are very attractive as they naturally provide 
a description of the dynamics at the population level and the resulting equations can be 
analysed using a wide range of analytical and numerical techniques. 

Particle-based simulations resolve the dynamics of each individual swimmer and from 
their positions and orientations, construct a picture of the dynamics of the suspension as 
a whole. As discussed in [T6|, particle-based models provide opportunities to (i) test con¬ 
tinuum theories, (ii) analyse hnite-size ehects resulting from a discrete number of swim¬ 
mers, (iii) explore more complex interactions between swimmers and/or boundaries, and in 
some cases, (iv) reveal the ehects of short-range hydrodynamic interactions and/or steric 
repulsion. Various models have been proposed in this context, each using diherent approx¬ 
imations to address the difficult problems of resolving the hydrodynamic interactions and 
incorporating the geometry of the swimmers. Some of the first such models used point 
force distributions to create dumbbell-shaped swimmers [211 1221 I2E], slender-body theory 
to model a slip velocity along the surfaces of rod-like swimmers [21112S|, or the squirmer 
model [2E1122] to examine the interactions between spherical swimmers [2E]- These initial 
studies provided important fundamental results connecting the properties of the individual 
swimmers to the emergence of collective dynamics. Based on their success, these models 
have been more recently incorporated into a number of numerical approaches for suspension 
and fluid-structure interaction simulations including Stokesian dynamics [221 EDI El], the 
immersed boundary method [321 ED], Lattice Boltzmann methods [SIES], and hybrid finite 
element/penalization schemes [36]. This has allowed for both increased swimmer numbers 
as well as the incorporation of other effects such as steric interactions, external boundaries, 
and aligning torques. 

In this paper, we introduce an extension of the force coupling method (FCM) [271 EH], 
an approach for the large-scale simulation of passive particles, to capture the many-body 
interactions between active particles. FCM relies on a regularized, rather than a singular, 
multipole expansion to account for the hydrodynamic interactions between the particles. 
It includes a higher-order correction due to particle rigidity by enforcing the constraint of 
zero-averaged strain rate in the vicinity of each particle. Since the force distributions have 
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been regularized, the total particle force, including that associated with the constraints, can 
be projected onto a grid over which the fluid flow can be found numerically. This allows the 
hydrodynamic interactions for all particles to be resolved simultaneously. This mesh can be 
structured and simple such that an efficient parallel Stokes solver can be used to hnd the 
hydrodynamic interactions. 

We extend FCM to active particles by introducing the regularized singularities in the 
FCM multipole expansion that have a direct correspondence to the surface velocity modes 
of the squirmer model |Z7|. With these terms included, we then rely on the usual FCM 
framework to resolve the hydrodynamic interactions in a very efficient manner. We show 
that by using the full capacity afforded by FCM, we are able to accurately simulate ac¬ 
tive particle suspensions in the semi-dilute limit with 0(10^ — 10®) swimmers. Using this 
method, we examine the influence of domain size on the steady-state polar order observed 
for squirmer suspensions. At the same time, we show that our method is quite versatile, 
being able to handle time-dependent swimming gaits, ellipsoidal swimmer shapes, and steric 
interactions, each at a minimal additional computational cost. We explore in detail how to 
incorporate biologically-relevant, time-dependent swimming gaits by tuning our model to the 
recent measurements of the oscillatory flow around Chlamydomonas Rheinardtii [SH] . These 
experiments revealed that considering time-averaged flows for such micro-organisms may 
oversimplify the hydrodynamic interactions between neighbors. Time-dependency is also 
closely associated with the way zooplankton feed, mix the surrounding fluid, and interact 
with each other Pi- As stated in [10], modelling micro-swimmers with a time-dependent 
swimming gait might be more realistic and should be included in mathematical models 
and computer simulations. We show that time-dependence can indeed affect the overall 
organization of the suspension. 

We organize our paper as follows: In section]^ we review FCM and present the theoretical 
background for its adaptation to active particles. Section details the numerical method, 
its algorithmic implementation and how the computational work scales with the particle 
number. In Section we validate the method and test its accuracy by comparing flow 
helds, trajectories, and pairwise interactions with previous results available in the literature. 
The effectiveness of our approach is demonstrated in Section where we present results 
from large-scale simulations of active particle suspensions. Finally, extensions of FCM to 
more complex scenarios are introduced in Section p We simulate suspensions of spheroidal 
swimmers and demonstrate the new implementation of time-dependent swimming gaits. 
Here, we also present preliminary results showing the effect of time-dependence on suspension 
properties. 

2. Squirmers using FCM 

The force-coupling method (FCM) developed by Maxey and collaborators |37l EH] is an 
effective approach for the large-scale simulation of particulate suspensions, especially for 
moderately concentrated suspensions at low Reynolds number. In this context, it has been 
used to address a variety of problems in microfluidics HU, biofluid dynamics m, and micron- 
scale locomotion ^31 ITTl H5] . FCM has also been extended to incorporate hnite Reynolds 
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number effects [10], thermal fluctuations m, near contact lubrication hydrodynamics [IHl 
Sn], and ellipsoidal particle shapes |SU]. With these additional features, FCM has been used 
to address questions in fundamental fluid dynamics in regimes where inertial effects are 
important and/or there is a high volume fraction of particles |l9l|51]. At the the same time, 
FCM has been used to address problems of technological importance, such as micro-bubble 
drag reduction |1B] and the dynamics of colloidal particles m- In this section, we expand on 
[52] and develop the theoretical underpinnings of FCM’s further extension to active particle 
suspensions using the squirmer model proposed by Lighthill [26], advanced by Blake iza. 
and employed by Ishikawa et al. [28|. To begin this presentation, we give an overview of 
FCM, establishing also the notation that will be used throughout the paper. 


2.1. FCM for passive particles 

Consider a suspension of Np rigid spherical particles, each having radius a. Each particle 
n, (n = 1,..., Np), is centered at Y„ and subject to force and torque Tn- To determine 
their motion through the surrounding fluid, we hrst represent each particle by a low order, 
hnite-force multipole expansion in the Stokes equations 

Vp-r^V^U = y^F„A„(x) ^Tn X V0n(x) Sn • V0n(x) 

n 

V-u = 0. (1) 


In Eq. Q, Sn are the particle stresslets determined through a constraint on the local 
rate-of-strain as described below. Also in Eq. ([^ are the two Gaussian envelopes, 

A„(x) = 

0„(x) = (2) 

used to project the particle forces onto the fluid. 

After solving Eq. (|^, the velocity, V„, angular velocity, 17^, and local rate-of-strain, E„, 
of each particle n are found by volume averaging of the resulting fluid flow. 



(3) 

(4) 

(5) 


where the integration is performed over the entire domain. In order for Eqs. (§-(§ to 
recover the correct mobility relations for a single, isolated sphere, namely that V = F/ (Gvrap) 
and Cl = T / {Sira^T]), the envelope length scales need to be a a = aj and gq = aj 
As the particles are rigid, the stresslets are found by enforcing the constraint that = 0 
for each particle n [38] . 
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2.2. Squirmer model 

In addition to undergoing rigid body motion in the absence of applied forces or torques, 
active and self-propelled particles are also characterized by the flows they generate. To 
model such particles, we will need to incorporate these flows into FCM. We accomplish this 
by adapting the axisymmetric squirmer model [2^ YFA |2S] to the FCM framework, though 
we note that a more general squirmer model [53] with non-axysimmetric surface motion 
could also be considered. 

The squirmer model consists of a spherically shaped, self-propelled particle that utilizes 
axisymmetric surface distortions to move through fluid with speed U in the direction p. If 
the amplitude of the distortions is small compared to the radius, a, of the squirmer, their 
effect can be represented by the surface velocity, v(r = a) = Vri -f vgO where 


OO 



( 6 ) 


OO 



(7) 


n=l 


Here, Pn{x) are the Legendre polynomials. 


VnicosO) = — -- sin 6*P'(cos 6*), 

n[n + 1) 


( 8 ) 


the angle 6 is measured with respect to the swimming direction p, and P^ix) = dPn/dx. In 
order for the squirmer to be force-free, we have 




(9) 


Following Ishikawa [28], we consider a reduced squirmer model where = 0 for all n and 
= 0 for all n > 2. We therefore only have the hrst two terms of the series. For this case, 
the resulting flow held in the frame moving with the swimmer is given by 




( 10 ) 


a 


4 


,4 


B2V2{cos9) 


( 11 ) 


r 
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where we have used U = 2i?i/3 from Eq. In terms of p, x, and r, this becomes 
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(12) 


= USj + US2 


Fig. shows an example of a flow held given by Eq. ( [I^ , as well as the hows uand u B 2 
related to the Bi and B 2 contributions. 


U 




Figure 1: Decomposition of squirmer velocity held for (3 = 1. 


While we have already seen that Bi is related to the swimming speed, it can be shown 
28] that B 2 is directly related to the stresslet 


G = (3pp - I) B2 


(13) 


generated by the surface distortions. This term sets the leading-order how held that decays 
like r“^. We can introduce the parameter (3 = B 2 /B 1 which describes the relative stresslet 
strength. In addition, if /3 > 0, the squirmer behaves like a ‘puller,’ bringing huid in along p 
and expelling it laterally, whereas if /3 < 0, the squirmer is a ‘pusher’, expelling huid along 
p and bringing it in laterally. 

To adapt this model to the FCM framework, we hrst recognize that the how given by 


Eq. (12) can be represented by the following singularity system in the Stokes equations 


Vp - = G ■ V 5(x) + — V^^(x) + HV^5(x 


6 


V-u = 0 


(14) 

(15) 
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where the degenerate quadrupole is related to Bi through 


H = —7r?7a^i?ip. 
3 


(16) 


and the stresslet G is given by Eq. (13). We can draw a parallel between these singularities 


and the regularized singularities used with FCM. The stresslet term in Eq. (14) is the 


gradient of the singularity system for a single sphere subject to an applied force. Accordingly, 
the corresponding regularized singularity in FCM is VA(x), where A(x) is given by Eq. (|^. 
It is important to note that even though we replace two singular force distributions with 
the one regularized FCM distribution, the particular choice of A(x) will yield flows that are 
asymptotic to both singular flow helds [2Z]. For the degenerate quadrupole, however, there is 
not a corresponding natural choice for the regularized distribution. Following [5l], we choose 
a Gaussian envelope with a length-scale small enough to yield an accurate representation 
of the singular flow, but not so small as to signihcantly increase the resolution needed in 
a numerical simulation (see Section 4.1). We therefore employ the FCM envelope for the 
force dipole and replace the singular distribution by V^0(x). Thus, for a single squirmer, 
the Stokes equations with the FCM squirmer force distribution are 


= G ■ VA(x) + HV^0(x) 

V-u = 0. 


(17) 


As we show Section]^ the flow satisfying Eq. (17) closely matches that obtained using the 
original singularity distribution, Eq. (14). 


2.3. Squirmer interactions and motion 

Using FCM, the task of computing the interactions between squirmers is relatively 
straightforward. We now consider Np independent squirmers where each squirmer has swim¬ 
ming dipole Gn and degenerate quadrupole H„. The squirmers may also be subject to 
external forces F„ and torques Tn- Using the linearity of the Stokes equations, we combine 
Eqs. ([^ and 0 to obtain 


Vp-r^V^u = ^F„A„(x) -h-T„ X V0„(x)-h Sn ■ V0n(x) 


+G„ ■ VA„(x) + H,V20 „(x) 

V-u = 0 


(18) 

(19) 


for the flow held generated by the suspension. 

After Ending the how held, we determine the motion of the squirmers using Eqs. ([^ - 
(|^ with two modihcations. First, we need to add the swimming velocity, t/p„, to Eq. 
Second, we must subtract the artihcial, self-induced velocity and the local rate-of-strain due 
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to the squirming modes. The self-induced velocity is given by 


W, = / A ■ H„A(x)d3x 


where the tensor A, in index notation, is m 


Aijix.) = 


Airrjr^ 


Sii 


SXiXn 


1 

V 




X^Xj 

- V 


erf 


(.Tev^) 




e(x) 


The self-induced rate-of-strain is given by 

K„ = / i (VR. G„ + (VR. G„)’') e(x)ci’x 


(20) 


( 21 ) 

( 22 ) 


(23) 


where the expression for the third rank tensor R can be found in [3H]- Taking these self- 
induced effects into account, the motion of a squirmer n is given by 

= Upn-Wn + j uA„(x)(i='x (24) 

fin = ^ y [V X u] 0„(x)d^x (25) 

= -K„ + iy [Vu+(Vu)^]0„(x)d3x. (26) 


As before, the stresslets S„ due to squirmer rigidity are obtained from the usual constraint 
on the local rate-of-strain, namely = 0 for all n. Squirmer positions Y„ and orientations 
Pn are then updated with the Lagrangian equations 


dYn 

dt 


V. 


(27) 


= fin X Pn. (28) 

In Section]^ we show through a comparison with the boundary element simulations from 
[28] that our FCM squirmer model recovers the velocities, angular velocities, stresslets (S^) 
and trajectories for two interacting squirmers for a wide range of separations. 


2.3.1. Including additional features 

With FCM, additional effects can readily be incorporated into the squirmer model and in 
our subsequent simulations, we consider several of them to demonstrate the versatility of our 
approach. For example, forces due to steric interactions and external torques experienced 
by magnetotactic or gyrotactic organisms can be considered by including them in F„ and 
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in Eq. (18). Additionally, we can extend the FCM squirmer model to ellipsoidal shapes by 
using the ellipsoidal FCM Gaussian distributions. Such a model can be carefully tuned by 
comparing with results from |55] and [56] • The effects of particle aspect ratio on suspension 
properties can then be explored systematically while still accounting for particle size effects 
such as Jeffery orbits. Finally, we are not limited to constant values for Bi and B 2 . By 
allowing these parameters to be functions of time, the FCM squirmer model can be used to 
explore how the swimmers’ strokes affect overall suspension dynamics. These extensions are 
addressed in Section [6l 


3. Numerical Methods 


3.1. Fluid solver 

The smoothness of the Gaussian force distributions allows FCM to be used with a variety 
of numerical methods to discretize the Stokes equations. It has been implemented with 
spectral and spectral element methods |l2l [501119| and finite volume methods mm in 
both simple and complex domain geometries. 

Here, we use a Fourier spectral method with Fast Fourier Transforms (FFTs) to solve 
the Stokes equations, Eqs. (18) - (19), for the fluid flow in a three-dimensional periodic 
domain. We set the zeroth Fourier mode of the velocity to zero, u(k = 0) = 0, to ensure no 
mean fluid flow through the unit cell. As the fluid encompasses the entire domain, this is 
equivalent to adding a uniform pressure gradient to balance the net force the particles exert 
on the fluid in order to maintain the system in equilibrium [521 EZj. This also ensures we have 
a convergent evaluation of the stresslet interactions [60] . As the squirmers can move without 
exerting a force on the fluid, we do permit a net material flux of squirmers through the unit 
cell. This corresponds directly to the far-held Stokesian Dynamics computations performed 
in [20] for squirmer suspensions. The FFTs are parallelized with the MPI library P3DFFT. 
This library uses 2D decomposition of the 3D domain, introducing a better scalabilty than 
FFT libraries that implement a ID decomposition. This decomposition has shown good 
scalability up to Nc = 32, 768 cores in Direct Numerical Simulation (DNS) of turbulence 
EH. 


3 . 2 . Computational work 

As explained in [19], the number of hoating point operations for FGM scales linearly 
with the number of particles, Np. We hnd the same scaling for our implementation of FGM. 
Fig. 1^ shows the computational time per times-step for Np up to 80, 000 particles with 
384^ ~ 6 ■ 10 ’^ grid points and Nc = 256 cores. 


3.3. Steric interactions 

Including steric repulsion is straightforward with FGM. These forces are introduced to 
both prevent particles from overlapping during the hnite time-step and to account for contact 
forces. For spherical particles, we use the steric barrier described in [US]- For particles n 
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4^v 



Figure 2: Scaling of the FCM with the number of squirmers Np. Computational time per 
time-step versus Np^ or equivalently, versus volumetric fraction Nc = 256 cores work in 
parallel for a cubic domain with 384^ grid points. 

and m, let — Y„ and r^m = llrnmll- The repulsive force experienced by n due to 

steric interactions is 



Fref 

' r 2 _ ' 

^ ref ' 

27 

^nrri') for '^ref’ 


F^ = 

2a 

R^ r - 4a2 
ref 

( 29 ) 


^0, otherwise. 




where is the magnitude of the force, the cut-off distance sets the distance over 
which the force acts, and the exponent 7 can be adjusted to control the stiffness of the force. 
Unless specihed, all the simulations are run with Qi^rjaU = 4, = 2.2a and 7 = 2 . 

From Newton’s third law, we obtain = —F'^. It is important to note that this force 
barrier is generic and it is not intended to model any particular physics. While we could 
instead use an exponential DLVO-like potential [HE], the force barrier can be evaluated more 
rapidly than an exponential potential and by adjusting its parameters we can incorporate 
repulsion of a similar strength and length scale as DLVO forces. A thorough study on the 
effect of this force barrier on the dynamics of particulate suspensions is provided in [U^ . 
Steric repulsion for spheroidal particles requires different modelling. In our study, steric 
forces and torques are introduced by using a soft repulsive potential with the surface-to- 
surface distance approximated by the Berne-Pechukas range parameter [63] . 

Using a direct pairwise calculation, the evaluation of steric interactions between all par¬ 
ticle pairs at each time step would require 0{Np/{2Nc)) computations per core. This cost 
is much greater than the 0{Np) cost of the hydrodynamic aspects of FCM. Therefore, 
instead of a direct calculation, we use the linked-list algorithm described in [ 6 l]- This 
method divides the computational domain into smaller sub-domains into which the parti- 
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cles are sorted. The edge-length of each sub-domain is slightly larger than These 

sub-domains are distributed over the cores where the steric interactions are evaluated. For 
a homogeneous suspension, this results in a per core cost for steric interactions that is 
O (l4iV2/(2iV,iV,)), where iV, is the total number of sub-domains. Since for large system 
sizes, we have Ng ~ ^ 14 and A'^c S> 1, the linked-list algorithm for steric interac¬ 

tions is much more efficient than the direct computation. 


3.4- Algorithm 

We summarize the overall procedure to simulate large populations of microswimmers in 
Stokes flow with the FCM: 


Initialize particle positions and orientations p 


JO) 


Start time loop: for A: = 1,..., Wt 

1 . Compute Gaussians Ai^^(x) and 0n^^( 


2. Update swimming multipoles Eq. (|13|), and Hn'"'', Eq. (|16|), which both 


), Eq. (§, 


ik') 

depend on pk , 


r{^) 


3. Compute steric interactions, Eq. (29), with the linked-list algorithm. 


4. Add additional forcing if any (gyrotactic torques, magnetic dipoles,...). 


5. Project the Gaussian distributions onto the grid (RHS of Eq. (19)), 


6 . Solve Stokes equations, Eqs. (18) - (19), to obtain the fluid velocity field u*^^)(x). 


7. Compute particle rate of strains Eq. (|26|), 


8 . If ||En^^|| > £, compute stresslets following 


a) Project all the multipoles onto the grid (RHS of Eq. (19)), 


b) Solve for Stokes equations, Eqs. (18) - (19), to obtain the fluid velocity field 
u(^)(x), 

9. Compute particle velocities Yn \ Eq. (24), and rotations Eq. (25), 


obtain and pn 


(fc+i) 


10. Integrate Eqs. (27) and (28) using the fourth-order Adams-Bashforth scheme to 


4. Validation 

f.l. Comparison with Blake’s solution 


As explained in Section |2.3[ in FCM the velocity held around an isolated squirmer n in 
an inhnite, quiescent huid is given by 


u 


FCM 


= A + R : G, 


(30) 


where A is a second rank tensor given by Eq. 
found in |3H1. Fig. shows, Blake’s solution Eq. 


(21) and the third rank tensor R can be 


(12), and the velocity held provided by Eq. 


(30) with two diherent values for the degenerate quadrupole Gaussian envelope size, a© and 


cre/2, that appears in the tensor A. The streamlines are identical, except for a near-held 
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Figure 3: Velocity field n/U around a puller squirmer (/9 = 1) swimming to the right, a) 
Blake’s solution; b) FCM solution with cxe for the degenerate quadrupole envelope, c) FCM 
solution with as/2 for the degenerate quadrupole envelope. 


recirculating region that appears when the width of the degenerate quadrupole envelope is 
(Te (Fig. 3b). In this region, however, the magnitude of the velocity is small compared 


to the swimming speed. This region should not signihcantly impact the squirmer-squirmer 
hydrodynamic interactions that we are aiming to resolve. 

A quantitative comparison of the velocity held is provided in Fig. The agreement 
with Blake’s solution is very good for r/a > 1.25 when using cxe for the width of the 
degenerate quadrupole envelope. As shown in Fig. 4b, the smaller envelope size (cre/2) 


matches Blake’s solution more closely for r/a < 1.2, with clear improvement at the front 
and rear of the squirmer. For this envelope size, the velocity held induced by the degenerate 
quadrupole matches exactly the analytical solution in Eq. (12) (not shown here). 

Below this width no quantitative improvement is observed as the remaining error comes 
from the dipolar contribution, \ib 2 - 

While this quantitative comparison provides a nice way to choose the degenerate quadrupole 
envelope size, we must also keep in mind the computational cost associated with decreasing 
this length scale. Even though it would yield a how held slightly more in register with 
Blake’s solution, resolving the length scale (Te/2 in a 3D simulation would require a grid 
with 8 times as many points as that needed for ue, the smallest length-scale already in FCM. 
This would increase computation times by at least an order of magnitude. In addition, the 
FCM volume averaging used to determine the squirmer translational and angular velocities, 
Eq. ([^ - ([^, will reduce the contribution of the localized velocity held discrepancies to the 
squirmer-squirmer interactions. In our subsequent simulations, we therefore utilize (Tq for 
the degenerate quadrupole envelope size since reducing this length-scale would signihcantly 
increase the computational cost, but only provide a minimal improvement. 


4-2. Pairwise interactions of squirmers 

In |2H], the authors performed a variety of simulations using the Boundary Element 
Method (BEM) to compute to high accuracy the pairwise interactions between a squirmer 
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, Blake 


Figure 4: Comparison with Blake’s solution, Eq. (12). a) Normalized difference, ||u" 
uFCM||^^^ between FCM and Blake’s solution for a puller squirmer (/3 = 1). The half-width 

for the degenerate quadrupole envelope is cr© . - : 10% iso-value, b) Velocity profile 

along the swimmer axis for Blake’s solution and the FCM approximation for the two different 
degenerate quadrupole envelope sizes. 


and an inert sphere, and between two squirmers. Here, we consider the same scenarios as 
those authors and compare results from our FCM simulations with their BEM results. 


4-2.1. Interactions between a squirmer and an inert sphere 

We first consider the interactions between an inert sphere (labelled “2”) located at a 
point r from the center of a puller squirmer (labelled “1”) with /3 = 5. The direction r/r 
forms the angle 9 with the swimming direction p of the squirmer. The problem setup is 
depicted in Fig. 

Fig. § compares the velocity of the sphere obtained using FCM simulations with the 
BEM results and far-field analytical solutions from |28]. We computed the far-field velocity 
in Fig. 6a using the expressions provided in [28] since it was not plotted in Fig. 6a of |28| for 
rja< 2.5. We suspect that the BEM results were also not plotted for this range. As FCM 
also resolves the mutually induced particle stresslets, it provides a more accurate estimation 
than far-held approximation of [28]. As a result, we see that the FCM results very closely 
match BEM, even in the range where r/a < 3. Similar trends are observed for the angular 
velocity of the inert sphere (Fig. 6b) and the stresslet components (Fig. [^ 6d). These 


comparisons illustrate the accuracy of the results that can be obtained using FCM, which 
closely matches BEM but incurs a fraction of the computational cost. 
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Figure 5: Sketch showing the set-up of our computations of the interactions between squirmer 
“1” and inert sphere “2”. 

4-2.2. Trajectories of two interacting squirmers 

We compute the trajectories of two puller squirmers {(3 = 5) and compare the results 
with the BEM simulations of |2B]- One squirmer initially swims in the x-direction, pi = e^., 
and the other one in the opposite direction, p 2 = —e^. They are placed with initial sepa¬ 
ration distance 6y = la,..., 10a in the transverse direction and 6x = 10a in the x-direction. 
The problem set-up is depicted in Fig. Since the squirmers may collide, we also include 
steric interactions provided by the force barrier (29). 

As shown in Fig. the trajectories match very well with the BEM results for Sy > 
2a. When 6y = la, the collision barrier and near-held hydrodynamic interactions play an 
important role in determining the overall squirmer trajectories. Fig. shows the effect of 
the steric repulsion parameters and on the squirmer trajectories. The specihc 
values of these parameters are provided in Appendix B. We see that by varying the barrier 
parameters one can obtain trajectories that closely match the results of |28j . 

5. Simulation Results 

5.1. Large suspensions of swimming micro-organisms 

The squirmer model has been widely used to investigate both the behavior of single 
swimmers [65l |Ml EH EH] , as well as their collective dynamics and interactions |29l ISOl |Ml 
EH EH] • Simulations of suspensions revealed that the overall population dynamics depend 
strongly on the squirming parameter (3. In particular, when |/d| is small, the isotropic state 
for a periodic suspension has been shown to be unstable, and the suspension evolves to a 
polar steady-state with a non-zero value of the polar order parameter 



(31) 


14 







(c) (d) 


Figure 6: The a) Radial velocity jlFrgl, b) Angular velocity c) Stresslet component 

|5'a;a;,2|, and d) Stresslet component \Sxy, 2 \ for the inert sphere “2” at a distance r from a 
puller squirmer (/3 = 5). 
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Figure 7: Initial configuration of the squirmers in the trajectory simulations. 
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Figure 8: Trajectories of two squirmers swimming in opposite directions with transverse 
initial distance 6y = la, 10a. Lines: data from [2H]. Symbols: FCM results. 
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(a) (b) 

Figure 9: Influence of the collision barrier parameters on the squirmer trajectories for the 
case where initially 5y = la. a) Trajectories. Crosses: data from [28]. b) Force barrier 
prohles (see also Appendix B). The line styles correspond those showing the trajectories in 
a). - : barrier used in Fig. 


This instability has been studied numerically by [30] and [69] using Stokesian dynamics 
with Np = 64 swimmers for volume fractions 0^ = 0.01 —)■ 0.5. Using the Lattice-Boltzmann 
method, [Mj observed the same behaviour for Np = 2000 and 0^ = 0.1. 


5.1.1. Polar order parameter 

Using our FCM model, we study this instability and the resulting polar order of a 
squirmer suspension. In particular, we examine how the domain size affects both the growth 
rate of the instability and the hnal steady-state. The influence of domain size has not been 
addressed previously for squirmer suspensions though it has been observed in simulations 
of rod-like swimmers m- We perform simulations of semi-dilute suspensions (0^ = 0.1) 
of puller squirmers (0 = 1) in triply-periodic square domains with edge lengths ranging 
from L/a = 14 to L/a = 116. As the volume fraction 0^ is hxed, varying L/a increases 


of the number of swimmers from Np = 64 (as in 


) to Np = 37, 659. We initialize 


a homogeneous, isotropic suspension by distributing the swimmer positions uniformly in 
the domain and the swimming directions uniformly over the unit sphere. Depending on 
domain size, the simulations are run to hnal time tf = 1000 — 1500a/t/ with a time step 
At = 0.005a/U. Thus, each simulation requires between 2 x 10^ — 3 x 10® time-steps. 

Figure 10 illustrates the polar ordered state for a simulations with L/a = 116, Np = 
37,659. We quantify the degree of alignment of each swimmer n, p„ ■ (p), with the mean 
steady-state orientation (p) given by 


(P) 



(32) 
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Figure 10: Snapshots of the orientational state in a semi-dilute suspension (0^ = 0.1) 
containing Np = 37, 659 swimmers with /3 = 1. (p) is the mean steady-state orientation 
vector on the unit sphere dehned in Eq. (32). p„ • (p) quantihes the degree of alignment 
with the mean direction (p) of swimmer n. 



t = lOOOa/U 


where (p) lies on the unit sphere At t = 0, there is no clear mean orientation, while at 
t = lOOOa/tl, a signihcant proportion of particles are aligned with the mean direction (p). 
Fig. 11a shows P{t), Eq. (31), for the simulations with different domain sizes. For 


each domain size, we see that the suspension evolves from the initial isotropic state to one 
that has polar order. We observe, however, that the hnal value, P°°, of the polar order 
parameter depends on the domain size. We hnd that it decreases as L/a (and Np) increases. 
The data also shows that as L/a —)■ oo, P°° decays like {L/a)~^P (or Np^^"^) and reaches 
an asymptotic value of —)■ 0.452. These results would indicate that polar order should 

also arise in an unbounded suspension. We also observe that in the polar ordered state, the 
average swimming speed is 4% less than its value in isotropic state which itself is nearly 
the free swimming value. Indeed, we remind that in our simulations, a zero net flux of the 
suspension is imposed, therefore, when the swimmers are aligned, they are slowed down by 
the backflow. 

From our simulations, we may also analyse how L/a affects the time evolution of the 
instability. In Fig. |lla we clearly see that the time it takes to reach the hnal polar state 
increases with domain size. We analyze this data in more detail in Fig. now plotting 


it in semilogarithmic scale. For each case, we hnd that after an initial transient state, the 
instability grows exponentially and with a growth rate that is independent of the system 
size (see inset hgure). It would certainly be interesting to investigate if this result could be 
reproduced by a linear stability analysis of continuum model for a squirmer suspension. 
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(b) 


Figure 11: Polar order P{t) in a semi-dilute suspension (0^ = 0.1) of squirmer pullers 
{(3 = 1). a) Time evolution of polar order depending on the number of swimmers. — 0- : 
L/a = 14, Np = 64; ^ : L/a = 19, Np = 174; : L/a = 38, Np = 1, 395; : 

L/a = 58, Np = 4, 707; : L/a = 77, Np = 11,158; ^ : L/a = 116, Np = 37, 659. b) 

Steady state value depending on Np [L/aY- - : fit linear with 1/ (Inset): 

Dependence of with 1 / y/Wp. 



Figure 12: Characterization of the polar instability. (Main hgure): Time evolution of polar 
order in semilogarithmic scale suggests an exponential growth of the instability. (Inset): 
Growth rate of the instability for different Np (and L/a). 
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Figure 13: Orientational distribution on the unit sphere {4>,0). a) Distribution before the 
transition to polar order: tU/a = 90. b) Distribution once the polar ordered state is reached: 
tU/a = 1110. 


5.1.2. Orientational distribution 

We can examine the polar order in more detail by computing the orientational distri¬ 
bution \1'(6', 0, t) dehned over the unit sphere [^. Here, 9 = cos“^(p^) corresponds to the 
elevation angle while 0 = tajLi~^{py/px) gives the azimuthal angle. Fig. 13 shows \b(0,0,t) 
normalized by the isotropic distribution \bo = l/( 47 r) at times before and after the transi¬ 
tion to polar order for the case where Np = 11,158 and L = 78a. As expected, before the 
transition, the distribution is nearly uniform over the surface of the sphere. After the polar 
state is reached, we see that \b(6*, 0 ) is narrowly distributed around the mean direction (p). 
We note that the steady state mean direction depends both on the random initial seeding 
of swimmers and the domain geometry. We hnd that system hrst aligns along an arbitrary 
direction due to the random initial seeding of swimmers, but as time goes on, the mean 
director tends to align with the normal to one of the periodic boundaries. As a consequence, 
the steady state mean direction is biased by the domain shape that breaks radial symmetry. 
We would like to stress that, however, as demonstrated here and in [30], the polar ordering 
itself is a consequence of squirmer hydrodynamic and steric interactions rather than the 
domain shape. Fig. 


14 


shows the steady-state averaged distribution 4/(6*, 0)|(p) 


in the frame 

where the mean direHTon is given hj 6 = 0 and 0 = 0. We see that the resulting distribution 
is axisymmetric in that it does not depend on 0. This is not surprising as the flow held 
induced by a squirmer is axisymmetric. Orientational distributions could also be obtained 
from continuum models, as was done for swimmers in external how helds HOI. Our results 
could be compared with continuum models of squirmer suspensions, though, to the best of 
our knowledge, the continuum models that are currently available in the literature predict 
a stable isotropic state for spherical swimmer suspensions. 
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Figure 14: Time-averaged steady-state orientational distribution in the frame of the mean 
orientation vector 4/(6*, 0)|^p^. a) Distribution over the unit sphere, b) Q : distribution of 

elevation angle 6 averaged over azimuthal angle 0; - : uniform distribution To(6*) = l/vr; 

—B- : distribution of azimutal angle 0 averaged over elevation angle 6] - : uniform 

distribution To(0) = l/(27r). 


5.1.3. Spatial distribution 

Along with the orientational distribution, we also examine the spatial distribution of 
squirmers by computing the Voronoi tesselation for the set of points corresponding to the 
squirmers’ centers. We have used the C-I--I- library Voro-|--|- mi to perform the computation 
and determine the volume, Vy, of each Voronoi cell. Fig. 15 shows the evolution of the 


standard deviation, ay, of the Vy distribution for the case where Np = 11,158 and L = 78a. 
We see sudden increase in ay as the suspension transitions to polar order, corresponding to 
a widening of the distribution. The inset in Fig. 15 shows the time-averaged distribution 


before and after the transition. Before polar ordering occurs, the mode of the distribution 
is Vy ~ 34a^, which is slightly less than the average value (Vy) = Li^ /Np = 42a^ that one 
might expect. During the polar steady-state, we hnd that the mode decreases to Vy ~ 30a^, 
which is indicative of a slight clustering of the particles. 

To further examine the spatial distribution, we compute the steady-state pair distribution 
function g{r,6) which gives the probability of hnding a squirmer at a distance r = |r| and 
with elevation angle 6 = cos“^(r • p/r) from another squirmer that has swimming direction 
p. Fig. 16a shows g{r, 6 ) for the case where Np = 11,158 and L = 78a. We see clearly that 
g{r, 9) depends not only on r, but on 9 as well. At steady-state, for a given squirmer there is 
a signihcantly higher probability of hnding another squirmer in front of it (6* = 0) rather than 
behind it. If we integrate g{r, 9) over 6*, we obtain the radial distribution function shown in 
Fig. 16b We hnd that g{r) exhibits a peak at two radii from the swimmer surface. This 


peak suggests the existence of particle clusters whose sizes are greater than two individuals. 
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Figure 15: Distribution of Voronoi cell volumes. (Main figure): Time dependence of the 
standard deviation, ay, of the Voronoi volnme distribntion normalized by the mean valne 

{Vy)- (Inset): Time average of the Voronoi volnme distribntion before ( - ) and after 

( - ) the the polar order transition. 
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Figure 16: Steady state pair distribution fnnction a) g{r,9). b) g{r). 
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Figure 17: Correlations in squirmer orientation at steady state, a) Ip{r,9). b) /p(r) 


5.1.4- Orientational correlations 

The hydrodynamic and steric interactions between the squirmers also lead to correlations 
in their orientations. We compute the steady-state orientational correlation function inni 

/p(r) = (P(x) ■ P(x + r)), (33) 


where the brackets (■) denote the ensemble average that we compute by averaging over both 
time and squirmer pairs. Fig. |17a| shows the correlations in the frame of a squirmer located 


at the origin and with p = z and the dehnitions of r and 0 are identical to those in Section 
1.3, We hnd that the highest correlations occur near contact (r 2) at the angles 9 k, tt/S 
and 9 ^ Svr/d. Despite the overall polar order, we can identify two distinct regions around 
9 = 0 and 9 = tt where the orientations are uncorrelated. We do, however, find a positive 
value for the correlations, even far away from the origin. This is most evident in the 9- 
averaged correlation, Jp (r), shown in Fig. 17b which approaches a finite value Jp = 0.22 as 
r increases. This is consistent with the observed long-range polar order of the suspension. 
We also hnd that Jp (r) > 0 for all r. This again may be a result of the strong polar ordering 
of the suspension. 


6. Extensions to ellipsoidal swimmers and time-dependent swimming gaits 

In this section, we demonstrate the extension of our approach to both spheroidal swim¬ 
mer shapes and time-dependent swimming gaits. This illustrates the versatility of FCM 
while preserving good computational scalability and an accurate treatment of hydrodynamic 
interactions. 
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6.1. FCM for ellipsoidal particles 

To extend FCM to passive ellipsoidal particles, one simply modifies the Gaussian en¬ 
velopes in Eq. ([^ [50]. For example, taking the orthonormal vectors § 1 , § 2 , and §3 to be 
aligned with the ellipsoid semi-axes having lengths Oi, 02 , and 03 respectively, the Gaussian 
envelope corresponding to A„(x) is 


Af (x) = (27r) ^/^(crA;lO-A;20-A;3) 


aO-A-sj ^exp 


-(x- Yr 

2 V r 


’Q'^SaQ(x-Y„ 


(34) 


where (Taji = for i = 1,2,3, Q= (ei §2 63 ), and 


a 


Sa = 


-2 

A;1 

0 

0 


0 

.-2 
' A;2 


0 

0 

.-2 
' A;3 


(35) 


A similar expression is used for with cr©;* = a^/ {Q^/n) ' . Beyond this, the underlying 

algorithm of projecting the particle forces onto the fluid and volume averaging the resulting 
fluid flow remains unchanged. The constraint that E„ = 0 for each n is still used to find 
the stresslets. Thus, using FGM to compute the motion and hydrodynamic interactions of 
ellipsoids does not require any additional steps in the algorithm described in Section 3.4[ As 


explained in Section 3.3, the modelling of steric interactions differs from that for spherical 
particles and our simple force barrier could not be used. An extensive validation of FGM 
for ellipsoidal particles is presented in [50] where many of the classical results for ellipsoidal 
particles in Stokes flow, e.g. Jeffery’s orbits, are shown to be recovered exactly with FGM. 
We note that for very slender particles, FGM would not be appropriate and computations 
based on slender-body theories [72| should be performed. 


To extend FGM to active ellipsoidal particles requires adding the stresslet and possible 


potential dipole terms to the multipole expansion. As for spherical particles (see Section 2.3 
Eq. (20) and (23)), these additional multipoles will lead to artificial, self-induced velocities 


and local rates-of-strain. These effects must be subtracted away using the formula derived 
in Appendix A. It is worth reiterating that these are the only steps that need to be added 
to the FGM algorithm for passive particles to simulate active ones. 


6.2. Spheroidal swimmer simulations 

Previous particle-based simulations |2S1 E2] and continuum theory results [IHl UHl 120] 
predict an unstable isotropic state for suspensions of prolate spheroidal pushers. Using 
FGM, we simulate a dilute suspension, cfy = 0.05, of Np = 1500 spheroidal pushers with 
aspect ratio = 3. For these simulations, the stresslet parameter is i ?2 = —1.5 and the 

degenerate quadrupole is set to zero, Bi = 0. Steric interactions are included by using a soft 
repulsive potential with a ~ decay and surface-to-surface distance approximated by the 
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Figure 18: Simulation of a dilute suspensions, = 0.05, of prolate spheroidal pushers, B 2 = 
— 1.5, Bi = 0, with aspect ratio = 3. a) Snapshot at t = 3017/ai. Colors indicate 

the velocity magnitude normalized by the individual swimming speed, b) Time evolution of 
polar order P{t). -: isotropic value of the polar order parameter 1/= 0.026. 


Berne-Pechukas range parameter [63]. Fig. 18a shows a snapshot of the suspension where 
clusters of swimmers have velocities nearly 1.6 times larger than the isolated swimming 
speed. As in [251 Eg, we see that the isotropic state is not stable and the polar order 
parameter increases with time (Fig. 18b). We do not see as large an increase in P{t) as in 


due to the relatively low aspect ratio of our swimmers. The tumbling effect due to 
Jeffery orbits and the steric torques that tend to align adjacent swimmers are much lower 
than they are for rod-like swimmers. 


6.3. Suspension dynamics of time-dependent swimmers 

Here, we show how to incorporate time-dependence into our model. Using the procedure 
outlined in [73], we determine the time-dependent multipole coefficients Bi{t) and i? 2 (t) 
from the experimental data provided in [39]. We then perform suspension simulations that 
show time-dependence at the level of individual swimmers can affect the overall suspension 
properties. 

6.3.1. A single time-dependent swimmer 

Recent experiments [33] quantihed the periodic swimming gait and resulting flow held of 
the algae cell Chlamydomonas Rheinardtii. They extracted the swimming speed, the induced 
velocity held, and the power dissipation, showing also that all can be represented as peri¬ 
odic functions of time. A recent theoretical investigation [73] showed that these quantities 
could be reproduced using a multipole-based model. In their study, they considered three 
time-dependent multipoles: a stresslet, a degenerate quadrupole (or potential dipole) and a 
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“septlet.” The stresslet decays as whereas the degenerate quadrupole and the “septlet” 
decay like r~^. Here, we utilize only the stresslet and degenerate quadrupole terms and hnd 
that they are sufficient to reproduce the features measured by [39]. 

The measured swimming speed from [39], U(t), can be represented using the truncated 
Fourier series in time (see ra): 

U(t) = oo + Oi cos {cot) + 02 cos {2cot) + bi sin {cot) + 62 sin {2ut) (36) 


where u is the frequency of the swimming gait. The mean swimming speed over one beat 
period, T = 27i/co, is given by Oq = 49.54a.s“^, where a = 2.5/im is the radius of the 
microorganism. The values for the remaining coefficients are provided in the supplementary 
information of |S|. From Eq. we can immediately determine the time-dependent 

degenerate quadrupole strength 

B.(i) = ^U(t) (37) 

in order to preserve the instantaneous force-free condition. Unlike this term, there is more 
than one way to calibrate the stresslet strength, i? 2 (t). For example, one could determine 
B 2 {t) using the power dissipation measurements from [39] and Eq. ([^ from [71] 


nd(t) = gVrr/a {^Bi{tf + 4H2(f)^) 


(38) 


for the power dissipated by a squirmer. Using this approach, we found that our resulting 
flow held did not match the experimental results from [39]. We instead determine i? 2 (t) 
directly from the experimental how held by htting to the location of the moving stagnation 
point. This is the same approach adopted by [H]. We utilize three Fourier modes to describe 
the time evolution of B 2 


B2{t) = Co -F Cl cos (cut-t- (pci) + C2 cos (2ci;t-f (Pc2) + C3COS (3a;f-F (pcs) 
+Si sin {cot + Lfsi) + S 2 sin {2ut + (^^ 2 ) + S 3 sin {3ut + 99 ^ 3 ), 


(39) 


are htted manually. Appendix C 


where the amplitudes co,ci,si,... and phases (pci,^si, 
gives the resulting values of these parameters. 

The phase diagram in Fig. [19a 
found by [73]. We extract the average value of (3{t) = B 2 {t)/Bi{t) over one beat cycle 


shows the value of Bi versus B 2 and it is similar to that 


/ B^(t)IBi(t)dt = 0.1 


(40) 


which corresponds to a puller squirmer with a relatively small stresslet magnitude. Fig. 


19b shows the resulting power dissipation as determined from Eq. (38). It reaches a peak 


value at t/T 0.3, which coincides with the time at which the swimming speed reaches its 
maximum value. We note that unlike [39] our swimmers generate axisymmetric flow fields 
and we are considering a 3D periodic domain. Despite this, we achieve a qualitatively similar 
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Figure 19: a) (-Bi(t), B 2 (t)) phase diagram for one beat cycle, b) Power dissipation nd(t) 
over one beat cycle. - : FCM, Q • results from [39] . 


power dissipation profile with slightly greater values during the hrst half of the beat cycle. 

Fig. shows the flow field around our model of C. Rheinardtii at six different times 
during its beat cycle. These time points are chosen to correspond to those in Fig. 3 of [3^ . 
We achieve very similar streamlines and, by construction, the position of the stagnation 
point matches the experimental data very well. We also note that our flow field is similar to 
that given by the multipole model found in [73| even though we do not include the rapidly 
decaying “septlet” term in our model. These results illustrate that our properly tuned, 
time-dependent squirmer model can yield flow fields very similar to those of real organisms. 


6.3.2. Interactions between two model C. Rheinardtii 

We first consider pairwise interactions between two model C. Rheinardtii for two initial 
conhgurations, 5y = la and 2a (cf. Section 4.2.2). We introduce a phase shift, A(p, between 
their swimming cycles. When A(p = 0, the swimmers are synchronized, whereas for A 93 = tt 
they are completely out of phase. Fig. 21 shows the effect of this phase shift on the 
trajectories. For each case, we show only the trajectory of the swimmer labelled “2” in 
Fig. 0 We also provide the trajectories of steady swimmers with /3 = 0 and /3 = 0.1 
for comparison. Recall that fl = 0.1 corresponds to the average value of 82 ( 6 )/Bi{t) in 
our model. We see that the trajectories, including the final positions and orientations, do 
depend on A(p. We find, however, that these variations are small compared to the swimmer 
separation distance (see the inset of Fig. 21 for the relative trajectories in the frame of one 
swimmer). The trajectories are also close to those for steady squirmers with /3 = 0, 0.1. 
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Figure 20: Snapshots of the time-dependent flow held around a model Chlamydomonas. 
Background grey levels represent the natural logarithm of the norm of the velocity held in 
radii. ■ : position of the stagnation point given by the FCM model. ♦ : position of 
the stagnation point measured by |39] . (Insets) Swimming speed along the beat cycle. 
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Figure 21: Relative trajectories of two model C. Rheinardtii swimming in opposite directions 
with stroke phase shift A(p and initial transverse distance 5y = la and 2a. The grey levels 

lighten as A(p increases from Aip = 0 —)■ A(p = 77r/4. - : trajectory of two steady 

squirmers with (3 = D. - ; trajectory of two steady squirmers with [3 = 0.1. (Top) 

5y = la. (Bottom) 5y = 2a 
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Figure 22: Sketch of the simulations with time-dependent swimmers. Each swimmer n 
generates time-dependent flow disturbances with the same amplitude and frequency, but 
the phase may be different for each swimmer. 


6.3.3. Collective dynamics of time-dependent swimmers 

In this section, we present results from suspension simulations using our unsteady model 
for C. Rheinardtii. To the best of our knowledge, similar results have not appeared in 
the literature. We aim to show in this initial study that time-dependence of individual 
swimmers can affect their overall distribution. Specifically, we show that the distribution of 
beat phases affects the steady-state polar order of the suspension. Fig. shows a sketch 
of the simulations with time-dependent swimmers. 

We consider a suspension of Np = 1395 time-dependent swimmers distributed uniformly 
in a triply periodic domain. The volume fraction is 0^ = 0.1. We take the initial swimming 
directions to be distributed uniformly over the unit-sphere. For swimmer n, its gait is 
characterized by Bi{t + Aipn) and B 2 (t A(p„), where A(p„ is the value of its beat phase. 
We consider two different distributions of the beat phases. The phases are either uniformly 
distributed with Aipn ^ [0; 27r], or Aipn = 0 for all n such that all swimmers are synchronized. 
We run the simulations for 3700 dimensionless time units, corresponding to 4000 beat cycles. 

the mean velocity of C. Rheinardtii isU = 49.54a.s“^. 


As mentioned in Section 6.3.1 


Since 

we take the time-step At = 0.0025a/f/, the total time for our simulations correspond to 
1.5 X 10®At. 

Fig. 


23a shows the evolution of the polar order parameter for these two distributions of 


beat phase. Also shown is the polar order parameter for suspensions of steady swimmers with 
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(a) (b) 

Figure 23: Polar order parameter and orientational distribution for suspensions of steady 

and time-dependent swimmers. - : time-dependent swimmers with random phase;.: 

synchronized time-dependent swimmers;-: steady swimmers with /d = 0;- : steady 

swimmers with f3 = 0.1. a) Time evolution of polar order, b) Steady state orientation 
distribution around the mean director - : uniform distribution To(6 ') = I/tt; 


swimming speeds equal to the mean swimming speed for the time-dependent case and either 
/3 = 0 or /3 = 0.1. We hnd that when there is a uniform distribution of the beat phase, we 
achieve results that match the steady case with /? = 0.1. On the other hand, if the swimmers 
are synchronized, we hnd that the polar order parameter matches that for the steady case 
with = 0. Fig. 23b shows the steady-state orientational distributions about the 

mean director. Again, we see that the distribution for case of a uniform distribution of beat 
phases matches that for the steady case with (3 = 0.1, while the synchronized suspension 
has the same distribution as the /d = 0 case. These new results show that the distribution 
of beat phases can strongly affect the orientational order of the suspension. 


7. Conclusions 

In this study, we presented an extension of FCM to compute the hydrodynamic inter¬ 
actions between a large number (10"^ — 10®) of active self-propelled particles in semi-dilute 
suspensions. Our approach builds from the spherical squirmer model developed by Blake 
izg by including in the usual FCM force distribution the regularized singularities that cor¬ 
respond to the surface squirming modes. We have shown that our model readily allows for 
different swimming gaits, as well as spheroidal swimmer shapes, and can also account for 
effects such as steric repulsion. We demonstrated the accuracy of our model by comparing 
velocity helds, pairwise interactions, and trajectories with the analytical or numerical results 
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given by the full squirmer model. In our squirmer suspension simulations, we recovered the 
polar order found in other studies, but also quantihed the effects of domain size on the hnal 
steady-state. Specifically, we showed that for puller suspensions with /3 = 1, the final value 
of the polar order parameter decreases and the time needed to reach steady state increases 
as the system gets larger. We did, however, find that these quantities appeared to converge 
to an asymptotic value and that polar order should be present even for an infinitely large 
domain. The scale of our simulations allowed us to compile robust statistics and examine in 
detail the orientational distribution and local micro-structure of the suspension. By extend¬ 
ing the model to time-dependent swimming gaits, we also illustrated for the hrst time, to 
best of our knowledge, that time-dependence at the level of individual swimmers can affect 
the hnal steady-state distribution of the suspension. 

We see that the squirmer model within the FCM framework provides an effective compu¬ 
tational approach to simulate active suspensions, allowing for particle numbers that begin to 
approach continuum levels. As a result, we feel that our computational model can both be 
used to complement experimental research, as well as inform the development of continuum 
models. For example, we have shown that the time-dependency of swimming gaits can be 
readily included in our model by tuning it to available experimental data [SH]. By tuning 
our model to similar experimental data, but for a wider zoology of microorganisms [9], we 
could assess the possible differences in collective dynamics exhibited by different species, or 
even look into how one species might interact with another. In addition, our model could be 
an effective approach to address questions regarding the mixing of background scalar helds 
[33] or passive tracers dynamics [731 Ei in active, time-dependent suspensions. The effects 
of tracer Brownian motion mi can be included with squirmer model and the competition 
between tracer diffusion and advection due to the flows produced by the swimmers can be 
assessed. 
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Appendix A: Self-induced effects for spheroidal swimmers 

In this Appendix, we show how to compute the artificial self-induced effects due to the 
squirming modes (see Section 2.3) for the case of spheroidal swimmers described in Section 
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I6.1[ 


In the following calculations, we consider an ellipsoidal squirmer located at the origin 
Y = 0 with swimming direction p = ei. 

Self-induced velocity W due to the degenerate quadrupole H 

Using the FCM envelopes for ellipsoidal particles, the force-distribution corresponding 
to the degenerate quadrupole generated by the squirmer will be given by 

= HV20"“(x). (A.l) 

where we also have cr© = (a©;i, cre; 2 , cre; 3 ), the width of the Gaussian envelope 0®^*(x) for 
each of the semi-axes. The resulting fluid flow can be obtained in Fourier space as 

= (A.2) 

The self-induced velocity, W, will be given by volume averaging this fluid velocity against 
the monopole Gaussian envelope A®^^(x) for ellipsoidal particles. Accordingly, the expression 
for W is 

R3 

(A.3) 

R3 

where we write the width of the Gaussian envelope A®^^(x) for each semi-axes as cr^ = 
c’'A; 2 , c’'A; 3 )- This expression can be viewed as Wi = where 

§'"(k)A'“(k)<i“k. (A.4) 

R3 


is the FGM self-mobility matrix that relates the particle velocity to the degenerate quadrupole 
coefficient. For a spheroidal particle {o'a -,2 = c’'A ;3 and (j ©;2 = o'©; 3 ) and H in the direction 
ei, we need only to consider to obtain the self-induced effects. This mobility matrix 

entry can be written as 


= 


Svr^ 

I 

1 

Svr^ 


^ (^1 - 0"“(k)A"“(k)d3k, 


1 ^ . 


1 + x' 


(A.5) 


4;2 + kl + kl) 


d^k. 
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ai (j0;i (Ja;i . ,1 , . 1 i. i.- j ‘^0;i ‘^©;2 cre;3 

where 7 = — = - = - is the spheroid aspect ratio and x = “ “ 


7 r\ 1/3 

6 . 


0-2 O'e-2 O'A-2 


CA;! CrA;2 0'A;3 


If we introduce cylindrical coordinates {ki = z, k 2 = rcos6, fc 3 = rsin6*), Eq. (A.5) 
becomes 
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z rt dz r 1 /I A K 1 

Let t = , „ which gives z = . and — =-^ and Eq. (A.6) becomes 
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(A.7) 


After integration by parts with respect to r, we have 
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dt, (A-8) 


and eventually 




( 1 -/ 2 ) 


2^2 
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(A.9) 


From the integrals in the Appendix in [SU], this maybe be expressed compactly as 

= -2C{Io-h). (A.IO) 


( 3 /2 \ — 

3271^/^Ha ^.2 (-|^) ) and Jq and A are coefficients whose 

provided in the Appendix of 


expressions are 
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The self-induced velocity W is then simply 


W = 


(ATI) 


and its artificial effect can be subtracted away from the total volume average velocity as 
done in Eq. (24). 


Self-induced rate of strain K due to the squirming dipole G 

In the frame of the swimmer, the non-zero entries of the self-induced rate of strain K' 
are given by 

K'li = 

K'22 = Tfgg.^22^22 + Tf|^Q.^23^33 + TfJ’c-.f^2i^ii! (AT2) 

K'33 = M ^ q ^ 232^22 + 

is the symmetric fourth-order FCM mobility tensor relating the particle rate-of-strain 
E to the swimming dipole coefficient G [50]. Its components read 


<^^212 = <^;f313 
^Eg!2222 = ^SG;3333 


^EG^323 


-2C (/ 
G(/i- 



1 “ -^ 2 ), 

I 2 ), 

1 3 

Iq —I\ — ~I2 

U T 2 i ^2 

Ir Ir 1 

4/0+2^1-4 



(A.13) 


where C is given above and I 2 is a coefficient whose expression is detailed in the Appendix 
of [^. To obtain the self-induced rate of strain in the lab frame, K, one just needs to apply 
a rotation operator onto K': 

K = QK'Q^, (A.14) 

where Q = ( ei 62 63 ) is the rotation matrix of the swimmer. 


Appendix B: Repulsive force parameters 

The table below lists the parameter values used in our squirmer trajectory calculations. 


-^rgf 

QnrjaU 

Exponent 7 

— 2.2 

4 

2 

2.4 

6 

1 

— 3 

3 

5 

2.04 

5 

5 


Table B.l: Parameters for the contact forces in Section 


4.2.2 


Fig. 


9b 


35 















Appendix C: Fitted values for C. Rheinardtii 


This appendix provides the values used to tune the time-dependent squirmer model 


presented in Section 6.3.1 


Co 

Cl 

C2 

cs 

Si 

S2 

S3 

4.5347 

64.053 

-84.7192 

-10.4545 

91.4529 

-92.6420 

-5.9849 


Table C.l: Magnitudes of the Fourier modes used to describe i? 2 (t) in radii.s ^ 






^S2 


1.7373 

3.5761 

-0.9154 

0.1666 

2.0054 

-1.7125 


Table C.2: Phases of the Fourier modes used to describe B 2 {t) in rad. 
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